In order to convert map coordinates to lat/lon data, we can use the ps2ll function like so:
# To call the function we pass in x and y values, and can also specify
# parameters such as earth_radius (km), eccentricity, and meridian.
x = -75
y = 150
[lat, lon] = PyPMT.ps2ll(x, y, earth_radius=6378388, eccentricity=0.08199, meridian=-100)
print(f"latitude = {lat}, longitude = {lon}")
latitude = -89.9985, longitude = -126.5651
# We can also pass in several x / y pairs at once.
# generate random numbers
x = np.random.uniform(low=-2000000, high=0, size=10000)
y = np.random.uniform(low=-1000000, high=0, size=10000)
# Test the function
lat, lon = PyPMT.ps2ll(x, y)
# Print the first 5 values of latitude and longitude
print("Latitude:", lat[:5])
print("Longitude:", lon[:5])
Latitude: [-84.1369 -72.4196 -84.8288 -81.6191 -87.5038] Longitude: [-157.2603 -101.1587 -99.2047 -166.8781 -114.5053]
If instead we wanted to convert longitude / latitude pairs to map coordinates for a polar stereographic system, we can use the ll2ps function like so:
# Pass in 2 latitude / longitude pairs.
lats = [-80, -65]
lons = [90, 120]
# We can specify the same paramters as ps2ll. Here we specify meridian.
x, y = PyPMT.ll2ps(lats, lons, meridian= -100)
print("x:", x)
print("y:", y)
x: [ -189134.02760898 -1773130.87598709] y: [-1072632.37224989 -2113135.09159733]
Other Mathematical Functions¶
PyPMT has a plethora of other functions for various purposes.
Suppose we want to estimate freeboard height above sea level, from ice basal elevation, assuming hydrostatic equilibrium. To do so, we can use the base_2_freeboard function:
# In this example, we alter the default snow density (rhos) and snow thickness (Ts).
height = PyPMT.base_2_freeboard(-4.5, rhos = 200, Ts = 0.75)
print(f"Freeboard height above sea level (m):", height)
Freeboard height above sea level (m): 1.22
Or if we wanted to estimate ice thickness from freeboard height (also assuming hydrostatic equilibrium), we can use the freeboard_2_thickness function:
PyPMT.freeboard_2_thickness(0.8)
7.47
Or suppose we have data taken along flight lines and we want to see where the two paths intersect. We can use path_crossing_ps71 to check if the paths intersect and, if they do, where the intersection lies.
# EXAMPLE 1
lat1 = [-71.0, -72.0]
lon1 = [0.0, 10.0]
lat2 = [-71.0, -72.0]
lon2 = [10.0, 0.0]
print(PyPMT.path_crossing_ps71(lat1, lon1, lat2, lon2))
(array([-71.5826]), array([5.]))
We can also find the distance between two or more points. Using path_dist, we calculate the distance (km) from Berlin to London, and from London to Paris. In this code, we set a reference point as London's coordinates.
lat = [52.5200, 51.5074, 48.8566] # Latitudes of Berlin, London, and Paris respectively
lon = [13.4050, -0.1278, 2.3522] # Longitudes of Berlin, London, and Paris respectively
ref_point = [51.5074, -0.1278] # London
units = 'kilometers'
# Now, the first distance (2nd number) is the distance from Berlin to London, and the second distance is the distance (third number) from London to Paris.
output = PyPMT.path_dist(lat, lon, units, ref_point = ref_point)
print(output)
[0, 934.5233499288536, 343.9231200909892]
Alternatively, we can find the cumulative distance along a path in polar stereographic coordinates by calling path_dist_ps.
# last parameter is the reference point
PyPMT.path_dist_ps([52.5300, 51.5074, 48.566], [13.4050, -0.1278, 2.3522], [51.5074, -0.1278])
array([-7777173.99984168, 0. , 2763411.11797939])
PyPMT.path_dist_ps([50, 51], [30, 31])
array([ 0. , 1016101.65930695])
Suppose we want to transforms georeferenced (zonal and meridional) vector components to cartesian (polar stereographic) coordinate components. To do so, we can use uv2vxvy like so:
# Polar Stereographic coordinates
x = -20000
y = 40000
# Velocity components in the Polar coordinate system
u = -30
v = 40
vx, vy = PyPMT.uv2vxvy(x, y, u, v)
print(vx, vy)
-44.72137860397854 22.360641666991803
Or conversely, suppose we want to transforms polar stereographic coordinates to georeferenced (zonal and meridional) coordinates. For this we can use vxvy2uv like so:
lat = np.array([-90, -80, -70])
lon = np.array([-180, -170, -160])
vx = np.array([1.5, 1.5, 1.5])
vy = np.array([-0.3, -0.3, -0.3])
u, v = PyPMT.vxvy2uv(lat, lon, vx, vy)
print(f"u = {u}")
print(f"v = {v}")
u = [-1.5 -1.52930608 -1.51214497] v = [ 0.3 0.03497006 -0.23112243]